BEA-explore-viz

library(tidyverse)
library(sp)
library(tmap)
bea <- readRDS("./data/bea.RDS")
counties <- readRDS("./data/county.RDS")

Cash receipts from marketings

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = mark_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = liv_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from marketings")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = mark_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = liv_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from marketings")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = mark_receipts), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,500000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Cash receipts from marketings") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_mark <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = mark_receipts[YEAR == max(YEAR)] - mark_receipts[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_mark <- sp::merge(counties, delta_mark, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_mark) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-15000,0,30000,60000,100000,200000,400000,800000,16000000,32000000,64000000), palette = "RdBu", border.col = "white",
              title = "Change in receipts from\nmarketings 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Cash receipts from livestock and products

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = liv_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = liv_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from livestock")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = liv_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = liv_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from livestock")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = liv_receipts), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,250000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Cash receipts from livestock") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_liv <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = liv_receipts[YEAR == max(YEAR)] - liv_receipts[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_liv <- sp::merge(counties, delta_liv, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_liv) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-60000,-30000,-15000,0,15000,30000,60000,120000,2500000), palette = c("#67001F","#D6604D", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC", "#053061"), border.col = "white",
              title = "Change in receipts from\nlivestock 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

Cash receipts from crops

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = crops_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = crops_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from crops")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = crops_receipts, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = crops_receipts, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Cash receipts from crops")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = crops_receipts), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,250000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Cash receipts from crops") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_crops <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = crops_receipts[YEAR == max(YEAR)] - crops_receipts[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_crops <- sp::merge(counties, delta_crops, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_crops) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-10000,0,10000,20000,40000,80000,160000,320000,640000,5000000), palette = "RdBu", border.col = "white",
              title = "Change in receipts from\ncrops 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Other income

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = oth_income, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = oth_income, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Other income")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = oth_income, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = oth_income, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Other income")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = oth_income), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,75000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Other income") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_othinc <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = oth_income[YEAR == max(YEAR)] - oth_income[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_othinc <- sp::merge(counties, delta_othinc, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_othinc) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-1000,0,5000,10000,15000,20000,40000,80000,200000), palette = "RdBu", border.col = "white",
              title = "Change in other income\n1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Warning: Values have found that are higher than the highest break

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Government payments

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = gov_pay, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = gov_pay, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Government payments")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = gov_pay, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = gov_pay, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Government payments")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = gov_pay), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,30000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Government payments") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_govpay <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = gov_pay[YEAR == max(YEAR)] - gov_pay[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_govpay <- sp::merge(counties, delta_govpay, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_govpay) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-5000,0,2500,5000,10000,15000,20000,40000,60000), palette = "RdBu", border.col = "white",
              title = "Change in government\npayments 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Warning: Values have found that are less than the lowest break

## Warning: Values have found that are higher than the highest break

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Misc. income

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = misc_income, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = misc_income, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Misc. Income")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = misc_income, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = misc_income, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Misc. Income")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = misc_income), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,40000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Misc. Income") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_miscinc <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = misc_income[YEAR == max(YEAR)] - misc_income[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_miscinc <- sp::merge(counties, delta_miscinc, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_miscinc) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-100,0,2500,5000,10000,15000,20000,40000,200000), palette = "RdBu", border.col = "white",
              title = "Change in misc.\nincome 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Production expenses

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = prod_exp, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = prod_exp, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Production Expenses")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = prod_exp, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = prod_exp, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Production Expenses")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = prod_exp), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,400000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Production Expenses") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_prodexp <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = prod_exp[YEAR == max(YEAR)] - prod_exp[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_prodexp <- sp::merge(counties, delta_prodexp, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_prodexp) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-10000,0,10000,20000,40000,80000,160000,320000,4500000), palette = "RdBu", border.col = "white",
              title = "Change in production\nexpenses 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Feed purchased

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = feed, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = feed, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Feed purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = feed, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = feed, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Feed purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = feed), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,50000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Feed purchased") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_feed <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = feed[YEAR == max(YEAR)] - feed[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_feed <- sp::merge(counties, delta_feed, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_feed) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-15000,-5000,0,5000,15000,30000,60000,1000000), palette = c("#B2182B","#F4A582","#D1E5F0","#92C5DE","#4393C3","#2166AC","#053061"), border.col = "white",
              title = "Change in feed\npurchased 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Warning: Values have found that are less than the lowest break

Livestock purchased

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = livestock, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = livestock, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Livestock purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = livestock, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = livestock, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Livestock purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = livestock), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,35000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Livestock purchased") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in receipts
delta_livp <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = livestock[YEAR == max(YEAR)] - livestock[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_livp <- sp::merge(counties, delta_livp, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_livp) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-40000,-20000,-10000,0,10000,20000,40000,80000,1000000), palette = c("#67001F","#D6604D", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC", "#053061"), border.col = "white",
              title = "Change in livestock\npurchased 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Warning: Values have found that are less than the lowest break

Seeds purchased

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = seeds, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = seeds, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Seeds purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
 # geom_point(aes(x = YEAR, y = seeds, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = seeds, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Seeds purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = seeds), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,35000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Seeds purchased") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in seeds purchased
delta_seeds <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = seeds[YEAR == max(YEAR)] - seeds[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_seeds <- sp::merge(counties, delta_seeds, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_seeds) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-1000,0,1000,2500,5000,10000,20000,215000), palette = "RdBu", border.col = "white",
              title = "Change in seeds\npurchased 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Fertilizer & lime purchased

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = fert_lime, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = fert_lime, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Fertilizer & lime purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
 # geom_point(aes(x = YEAR, y = fert_lime, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = fert_lime, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Fertilizer & lime purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = fert_lime), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,50000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Fertilizer & lime purchased") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in fertilizer purchased
delta_fert <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = fert_lime[YEAR == max(YEAR)] - fert_lime[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_fert <- sp::merge(counties, delta_fert, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_fert) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-1500,0,2500,5000,15000,30000,60000,120000,500000), palette = "RdBu", border.col = "white",
              title = "Change in seeds\npurchased 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Warning: Values have found that are less than the lowest break

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Petrol purchased

#geom_point
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_point(aes(x = YEAR, y = petrol, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = petrol, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Petrol purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

## Warning: Removed 440 rows containing missing values (geom_point).

#geom_smooth
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  #geom_point(aes(x = YEAR, y = petrol, color = Region), alpha = 0.2) +
  geom_smooth(aes(x = YEAR, y = petrol, group = Region, color = Region), alpha = 0.2) +
  theme_bw() +
  labs(title = "Petrol purchased")

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

## Warning: Removed 440 rows containing non-finite values (stat_smooth).

# boxplot
bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  ggplot(.) +
  geom_boxplot(aes(x = as.factor(YEAR), y = petrol), outlier.shape = NA) +
  coord_cartesian(ylim = c(0,20000)) +
   theme(axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  labs(title = "Petrol purchased") 

## Warning: Removed 440 rows containing non-finite values (stat_boxplot).

# change in petrol purchased 
delta_petrol <- bea %>%
  filter(GEOID %in% unique(counties$GEOID)) %>%
  group_by(GEOID) %>%
  summarize(DELTA = petrol[YEAR == max(YEAR)] - petrol[YEAR == min(YEAR)]) 

## `summarise()` regrouping output by 'GEOID' (override with `.groups` argument)

delta_petrol <- sp::merge(counties, delta_petrol, by = "GEOID",duplicateGeoms = TRUE, all = TRUE)
#delta_mark$LOSERS <- ifelse(delta_mark$DELTA < 0, "Decrease", "Increase") 

tm_shape(counties) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(delta_petrol) + 
    tm_polygons(col= "DELTA", style = "fixed", breaks = c(-500,0,1000,2500,5000,20000,50000,140000), palette = "RdBu", border.col = "white",
              title = "Change in petrol\npurchased 1969-2019") +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

## Variable(s) "DELTA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.